function[x_min f_min] = gridsearch_avg(candidate_indices, parameters,  Y, vdjQ, k, side, retailer,  Npoints, method)

    %method = 1 uses GMM
    %method =2 uses ML 
    
global transform 

m1 = Npoints; %mu
m2 = Npoints; %sigma
m3 = Npoints; %c


MLmu = -3;
MUmu = 1;       
MLsigma = .1;    
MUsigma = 5;   
MLc = -10;      
MUc = 0;      

G1 = linspace(MLmu,MUmu,m1);  %mu
G2 = linspace(MLsigma,MUsigma,m2); %sigma
G3 = linspace(MLc,MUc,m3);  %c

x1_min = MLmu;  %0 lower bound 
x2_min = MLsigma;  %0 lower bound 
x3_min = MLc;   %0 lower bound 

f_min = 1e+10;

if retailer==0     
for i1=1:m1
    for i2=1:m2
        for i3=1:m3
           
        candidate_params = [G1(i1), G2(i2), G3(i3)];
        
        if method==1
             f = gmm_fun_y_c_avg(candidate_params, candidate_indices, parameters, Y, k, side, retailer,1,0);
        else
             f =ml_fun(candidate_params, candidate_indices, parameters, Y, k, side, retailer);
        end 
            if f<= f_min
            f_min=f;
            x1_min = G1(i1);
            x2_min = G2(i2);
            x3_min = G3(i3);
            end
        end
    end
end 

if transform==1 
    x2_min=exp(x2_min);
end 

x_min = [x1_min, x2_min, x3_min];
 
elseif retailer==1
    
 for i1=1:m1
    for i2=1:m2
   
        candidate_params = [G1(i1), G2(i2)];
        
        if method==1
             f = gmm_fun_y_c_avg(candidate_params, candidate_indices, parameters, Y,k, side, retailer,1,0);
        else
             f =ml_fun(candidate_params, candidate_indices, parameters, Y, side, retailer);
        end 
            if f<= f_min
            f_min=f;
            x1_min = G1(i1);
            x2_min = G2(i2);
            
            end
    end
 end 

if transform==1 
    x2_min=exp(x2_min);
end 

x_min = [x1_min, x2_min];
 
end 
    

   